function shooting

%  Solves the BVP using linear shooting
%       y'' + p(x)y' + q(x)y= f(x)   for xL < x < xr  '

%  where
%      y(xl) = yL  and y(xR) = yR

%  p=0, q=-1, f=sin(2*pi*x)
%  xL=0, yL=0  and  xR=1, yR=0

% clear all previous variables and plots
clear *
clf

% set boundary conditions and parameters
	xL=0; yL=0;
	xR=1; yR=0;

% calculate and plot exact solution
xx=linspace(xL,xR,100);
exact=sin(2*pi*xx)/(1+4*pi*pi);
plot(xx,exact,'k')
hold on

% define title and axes used in plot
xlabel('x-axis','FontSize',14,'FontWeight','bold')
ylabel('Solution','FontSize',14,'FontWeight','bold')
title('BVP: Linear Shooting with RK4','FontSize',14,'FontWeight','bold')

% have MATLAB use certain plot options (all are optional)
box on
axis([0 1 -0.03 0.03]);
set(gca,'ytick',[-0.03 -0.02 -0.01 0 0.01 0.02 0.03]);
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 

% loop used to increase N value
for in=1:3

	% set number of points along the x-axis
	if in==1
		n=2
	elseif in==2
		n=4
	elseif in==3
		n=8
	end

	% generate the points along the x-axis, x(1)=xL and x(n+2)=xR
	x=linspace(xL,xR,n+2);
	h=x(2)-x(1);

	% calculate and then plot solution
	y0=[1 0];
	y1=rk4('de_0',x,y0,h,n+2);
	y0=[0 1];
	y2=rk4('de_0',x,y0,h,n+2);
	y0=[0 0];
	y3=rk4('de_f',x,y0,h,n+2);
	b=(yR-yL*y1(n+2)-y3(n+2))/y2(n+2);
	y=yL*y1+b*y2+y3;

	% plot the solution
	if in==1
		plot(x,y,'--sr','LineWidth',1,'MarkerSize',7)
		legend(' Exact',' N = 2',3);
		set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold');
		pause
	elseif in==2
		plot(x,y,'--ob','LineWidth',1,'MarkerSize',7)
		legend(' Exact',' N = 2',' N = 4',3);
		set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold');
		pause
	elseif in==3
		plot(x,y,'--*m','LineWidth',1,'MarkerSize',7)
		legend(' Exact',' N = 2',' N = 4',' N = 8',3);
	end	
	% Set legend font to 14/bold    
	set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold'); 
end
hold off

function g=q(x)
g=-1;

function g=p(x)
g=0;

function g=f(x)
g=-sin(2*pi*x);


% right-hand side of DE
function z=de_f(x,y)
z=[y(2) -p(x)*y(2)-q(x)*y(1)+f(x)];

% right-hand side of homogeneous DE
function z=de_0(x,y)
z=[y(2) -p(x)*y(2)-q(x)*y(1)];

% RK4 method
function ypoints=rk4(f,x,y0,h,n)
y=y0;
ypoints=y0(1);
for i=2:n
	k1=h*feval(f,x(i-1),y);
	k2=h*feval(f,x(i-1)+0.5*h,y+0.5*k1);
	k3=h*feval(f,x(i-1)+0.5*h,y+0.5*k2);
	k4=h*feval(f,x(i),y+k3);
	yy=y+(k1+2*k2+2*k3+k4)/6;
	ypoints=[ypoints, yy(1)];
	y=yy;
end;